Total energy calculation of high pressure selenium: The origin of incommensurate 
modulations in Se-IV and the instability of proposed Se-II 
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We present calculation of the high pressure crystal structures in selenium, including rational 
approximants to the recently reported incommensurate phases. We show how the incommensurate 
phases can be intuitively explained in terms of imaginary phonon frequencies arising from Kohn 
anomalies in the putative undistorted phase. We also find inconsistencies between the calculated 
and experimental Se-II phase - the calculations show it to be a metastable metal while the experiment 
finds a stable semiconductor. We propose that the experimentally reported structure is probably in 
error. 

PACS numbers: 61.50.Ks,62.50.+p 
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The group VI elements selenium and tellurium are 
semiconductors with a hexagonal crystalline structure 
at ambient pressure. Both elements undergo a series of 
phase transitions as pressure is increased, transforming 
to various complex structures denoted by roman numer- 
als I, II, III, IV. The generally accepted crystal struc- 
tures associated with these phases in selenium is hexag- 
onal Se-I; monoclinic Se-II from 14GPa; triclinic Se-III 
from 23GPa; orthorhombic Se-IV from 28GPa; and sub- 
sequently simple higher symmetry phases^]. Tellurium 
undergoes a similar set of phase transitions at lower pres- 
sures, without the equivalent of Se-II. 

In 2003 new work using combined powder and single- 
crystal x-ray diffraction[2} found that the Te-III phase 
had a more complex structure than had previously been 
thought, involving a modulated distortion of the unit cell. 
It was suggested that Se-IV had a similar structure ex- 
perimental confirmation of this is in progress[3j. 

Density functional theory 0, 0, using pseudopoten- 
tials, plane waves and local exchange correlation func- 
tional has been the method of choice for theoretical 
study of selenium. Within the local density approxi- 
mation(LDA) the p-bonded covalent chain structure and 
anomalous compression mechanism (one lattice parame- 
ter increases with pressure), was demonstrated Jfj. The 
anisotropy showed that the LDA introduced an error 
equivalent to a pressure, which can be improved by in- 
cluding gradient corrections (GGA)0. Since then band 
structures of the complex high pressure phases have been 
investigated using the experimental parameters without 
full relaxation of the crystal structure or incommensurate 
modulation^, and the liquid phase has been investigated 
up to high pressures 0. 

Although there are several papers on individual crys- 
tal structures, proposed transition mechanisms, and the 
simple high symmetry structures at pressures above Se- 
IVjlfJ, we are not aware of published calculated phase 
transition pressures in SeI-IV0|- This may be because 
of their structural complexity, because the enthalpy dif- 



ferences between the phases are extremely small, or, be- 
cause of imperfect agreement with the reported sequence 
of transitions. 

In this letter we report density functional calcula- 
tions of the first four phases of selenium, identify- 
ing the stable structures over the pressure range to 
40GPa. Subsequently we investigate the nature of the 
incommensurately-modulated fourth phase of Se, approx- 
imating the incommensurate structure and showing it to 
be an extreme case of a Kohn anomaly. 

We used the computer program VASP^^h treating ex- 
change correlation effect with a GGA which is essen- 
tial to obtaining accurate pressures for the low density 
structure^ Q- The wavefunctions are sampled using 
the method of Monkhorst and Pack, and the energies 
converged to lmeV with respect to k-point density and 
0.2meV with respect to plane wave cutoff. For Se I, II and 
IV a k-point mesh of 15 3 was sufficient while for Se-III we 
used 17 3 . Each structure was fully relaxed with respect 
to lattice parameters and internal coordinates accord- 
ing to ab initio calculated forces and stresses. The core 
region is described with ultra-soft psuedopotentials |l3| . 
For consistency, we performed calculations on reported 
structures of all phases I through IV. To calculate band 
structures a self-consistent calculation using the afore- 
mentioned meshes was first carried out, to obtain the 
ground state charge density, and then the single parti- 
cle states of the Kohn-Sham Hamiltonian were evaluated 
non-self consistently on a finer mesh. . 

Se-I is the ambient pressure phase of selenium. It has 
a hexagonal unit cell, with the atoms arranged in he- 
lical chains running parallel to the &-axis. Each unit 
cell contains three atoms, in positions (0,0, u) (u, 4,0) 
(-it, | - it) with u=0.289 at 20GPa. The equilibrium 
lattice parameters at this pressure are a=c=3.569A; 
b=5.135A. In accordance with previous work0 0, 0], 
we find that as the pressure is increased, the size of the 
b parameter at first increases, then starts to contract at 
lOGPa. However this change is a tenth smaller than the 
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size of the reduction in the other two directions. 

Se-II is reported to have a monoclinic structure, with 
six atoms per unit cell, at positions (0,0,0), (u, 0, v), 
(— u, 0, —v) and at (|, |,0) relative to these. At 22GPa, 
u=0.147 and v=0.339. The equilibrium lattice param- 
eters of the cell at this pressure are a=6.562Afe = 
2.665lc=6.326A/3=105.3° This is a layered structure, 
with each atom closely bonded to four other atoms. On 
the increase of pressure, the parameter a reduces much 
faster than b and c as the layers are being pushed to- 
gether. The angle j3 changes very little over the calcu- 
lated pressure range. 

FigQ] shows the density of electron states of Se-II at 
2GPa and 22GPa. The effect of pressure is primarily 
to broaden the distribution. Consistent with previous 
work, there are two sets of bands, corresponding to 4p 
electrons^^. There is no band gap, showing that at 
either pressure the material is a metal|l6j. 

Se-III is reported to be triclinic, with six atoms per 
unit cell. At 28GPa, the equilibrium lattice parameters 
of the cell are =3.460 Ab=12.240 Ac=2.624 Aa=91.085° 
/3=113.445° 7=89.608°. Se-III also has very anisotropic 
compression, the largest decrease in the cell is in the a 
direction, which decreases by 10% over a range of 20GPa, 
compared with 5% for c and only 1.5% for b. Angles a, (3, 
7 change little as the pressure is increased. At 24GPa, the 
calculated atomic positions are: (-0.003,-0.001, -0.024) 
(0.465. 0.167. 0.389) (0.030, 0.331, 0.064) and ( ±, §, \) 
relative to these. Band structure calculations show that 
Se-III is metallic. 

We modeled the Se-IV cell by relaxing ions and 
lattice parameters from the primitive cell reported 
experimentally^ without the incommensurate distortion 
(Fig|2J. At 35GPa, the equilibrium lattice parameters 
were a=3.383A, b=4.07lA, c=2.613Aand = 115.644° 
This compares with the experimental results: a=3.299A, 
b=3.996A, c=2.587A, and /3 = 113.105°. The calculated 
unit cell parameters are about 2% larger than the exper- 
imental value at the same pressure. This is typical for 
the GGA. 

It is impossible to do calculations on an incommen- 
surate cell using periodic boundary conditions, however 
one can use a rational approximant assuming that a 
distortion close to the observed q will also be stable rel- 
ative the undistorted structure. Here we approximated 
the distorted structure with tripled, quadrupled or hep- 
tupled cells corresponding to q=g, \ and |. In each case 
the calculation was started with frozen-in values of the 
incommensurate wave displacement vector taking Te-III 
data u x = 0.0249 and u z = 0.1026 as initial symmetry- 
breaking displacements. 

For q=i calculations were done at four different pres- 
sures: 25, 30, 35, and 40GPa. In all cases, the displace- 
ment remained and, the enthalpy of each atom in the 
distorted structure is 3meV smaller, with the distorted 
structure being slightly larger than the undistorted one. 
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TABLE I: Details of frozen phonons at q~j and q=| at var- 
ious pressures: amplitudes of modulations in the xz plane 
for q=T correspond to maximum displacement u, interme- 
diate displacement v and two atoms fixed at unmodulated 
positions. For q=T we have chosen the frozen phonon to be 
in phase with the unmodulated cell: within error the ratio 
u/v — \p2 implies that the distortion can be described as 
a single harmonic phonon. For q=f the displacement pat- 
tern is again harmonic, u x and u z represent the amplitude 
of a fitted sine wave, but the phase of the modulation is not 
locked. Enthalpy differences are of order 2-3 meV/atom (10- 
30meV/cell) from the unmodulated SelV calculated with the 
same size, kpoint sampling and supercell. The close match 
between modulated and unmodulated cells means that the 
enthalpy differences are significant. 

Table|l] shows a summary of results found. The energy 
differences to the undistorted phase are small, but lie 
within k-point and cutoff convergence. We tested this 
to greater accuracy which can be obtained by eliminat- 
ing sampling errors for the crucial energy differences by 
comparing cells with and without the distortion, using 
identical k-point and cutoffs [lif. 

The calculated distortions at q = \ and q = I are 
smaller than those found experimentally at 35GPa. This 
is to be expected since the calculated distortions are not 
the most favored one. Likewise, the very small enthalpy 
difference (-3meV/atom) is a lower bound on the sta- 
bility of the incommensurate wave. The displacement 
eigenvector (transverse, u z /v z — 0.22 ± 0.04 at 40GPa is 
in excellent agreement with the experimental value (0.26 
at 42GPaQ). 

For q=^ we carried out calculations at six different 
pressures: 10, 20, 25, 30, 35, and 40GPa, the last three 
when Se-IV was predicted earlier to be stable. At all six 
pressures, the structure relaxed back to the unperturbed 
structure, as did calculations with q = =. 

While DFT is clearly a sufficient theory to explain the 
incommensurate modulation, a more intuitive explana- 
tion exists: modulation arises from a Kohn anomaly (in- 
teraction between nested Fermi surface and a phonon) 
so large that the phonon energy becomes negative and 
the soft mode "freezes in". Experimentally |2j, the in- 
commensurate q- vector drops from 0.31-0.28 between 30 
and 50GPa, remaining at 0.28 until 70GPa|3j. A stan- 
dard Kohn anomaly in the free electron picture would 
occur at q — 2Uf, and soften (i.e. reduce the frequency 
of) the phonon at that wavevector. Here the softening 
of the mode is so great that it freezes in (i.e. the fre- 
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quency is imaginary) , lowering the symmetry of the unit 
cell and allowing coupling to the unit cell parameters. 
The stability arises because the phonon mode opens a 
pseudogap at the Fermi energy, which can been seen as a 
dip in the electronic density of states: FigO shows that 
the distortion has exactly this effect. We examined the 
band structure of the undistorted cell and find that the 
calculated Fermi surface cuts the band in the (0, k, 0) di- 
rection at values of k rising from 0.1 at 20GPa to 0.26 
at 40GPa, while the distorted structure has a pseudogap 
along this direction. Although this does not give us the 
exact value of q from the 2/cp of the undistorted phase, 
the agreement is sufficiently close to be confident of the 
mechanism. 

Fig0| shows the phase stability deduced from our cal- 
culations of enthalpy of each phase from 1 to 40GPa. 
The stable structure of selenium at each pressure is that 
with the lowest enthalpy, and although the differences in 
enthalpy are small throughout the whole pressure range 
they are above our convergence errors. 

The most notable feature is that our structure for Se-II, 
reported experimentally at MGPafiH El Elj . is not pre- 
dicted to be stable at any pressure. It is higher in energy 
by around 8meV/atom than Se-III, which is greater than 
convergence of the calculation. This may be because of 
an entropic effect: a symmetry breaking distortion may 
have lower energy with Se-II a dynamically stabilized 
phase with a transition below room temperature |22l l2.ll] 
(8meV corresponds to a temperature of 100K). However, 
this would not explain the band structure observed in the 
calculation - along with previous authors[iEi| we predict 
that Se II is metallic while experimentally it is a semi- 
conductor. This is a strong indication that the atomic 
positions used in the calculation (starting from the ex- 
perimental ones) are incorrect and we conclude that there 
may be an error in the experimentally reported crystal 
structure 

We predict a phase transition from hexagonal Se-I di- 
rect to monoclinic Se-III at 21.5GPa accompanied by a 
volume reduction of 5.5%; experimentally Se-III is ob- 
served at 23GPa0, HJ. The transition from Se-III to 
Se-IV is predicted at 27.5GPa, compared to the experi- 
mental value of 28GPa[l9l|2lJ,|25j. These excellent agree- 
ments give us confidence that our method correctly de- 
scribes high pressure selenium. 

To summarize, ab initio total energy and band struc- 
ture calculations on four different phases of selenium are 
reported. The generally accepted sequence of phase tran- 
sitions is not reproduced although unit cell and internal 
parameters are in good agreement with experiment. The 
calculated instability of the experimentally reported Se- 
ll structure, combined with our prediction that it should 
be metallic, leads us to suspect that this phase may have 
been mischaracterized, and there is very recent evidence 
for unexplained diffraction peaks in Se-Il|26j|. 

Body-centered monoclinic Se-IV is unstable against a 




-15 -10 -5 

Energy (eV) 



FIG. 1: Calculated density of states for Se II structure from a 
25 3 k-point grid at OGPa (dotted line) and 22GPa (solid line) 
with ions and unit cell fully relaxed within the spacegroup 
reported experimentally. The Fermi energy has been shifted 
to zero, although it lies at a minimum in the density of states, 
no band gap is present and this structure would result in a 
metallic material 16]. 
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FIG. 2: Four cells of the modulated body-centered monoclinic 
cell of Se-IV viewed down the a axis. The dashed line shows 
the modulation at the experimental value of g=0.28, the cal- 
culations were done on cells with this modulation adjusted to 
?=| and q = 4 corresponding to a modulation wavelength 
of 3 and 4 cell respectively. Atoms labelled with | lie at that 
fractional coordinate in the c direction, others lie in the plane. 
(Figure courtesy of M.I.McMahon) 

phonon modulation with q = j or q = S. Inspection of 
the band structure of the undistorted phase leads us to 
conclude that this instability arises from a Kohn anomaly 
in the phonon spectrum of sufficiently large magnitude 
that the phonon eigenenergy is negative (i.e. the mode 
is unstable). The wavevector of the nested Fermi sur- 
face is incommensurate with the crystal structure and 
varies with pressure. This gives a theoretical picture of 
the recent unexplained experimental observation of an 
incommensurate Se-IV. 
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FIG. 3: Calculated Kohn-Sham band structure and density 
of states for Se IV at 35GPa with (solid line) and without 
(dotted line) the modulation at q — i. The Fermi energy, 
Ef = 6A6eV in both cases, is shown by the vertical line. The 
Kohn anomaly arises because the distortion perturbs states 
at (0, |,0) which lie close to Ef, leading to the drop in the 
density of states at Ef 
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FIG. 4: Calculated enthalpy differences between Se II and 
other phases. The stable structure at is the one with the 
lowest enthalpy at each pressure. Circles denote Se I, squares 
Se III and diamonds unmodulated Se IV. Errors are of order 
the size of the symbols 
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